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Abstract 

The  first  step  of  fuel  cell  stack  characterization  consists  in  experimental  testing  of  its  global  electrochemical  performances  under  various  working 
conditions.  The  second  step  is  devoted  to  link  the  collected  data  to  local  species  concentrations,  membrane  hydration  and  resulting  current  density 
in  the  stack.  In  the  project  PICOS,  we  have  to  characterize  a  1.5  kW  12-cells’  stack  of  330  cm2.  The  stack  performances  are  measured  on  a  specific 
test  bench  and  analyzed  using  a  fuel  cell-stack  modeling  tool.  The  test  bench  allows  to  control  the  oxygen  (or  air)  and  hydrogen  flow  rates  and 
humidity  at  the  inlet  of  the  stack,  the  pressure  at  the  outlet,  and  the  delivered  current.  Each  cell  potential  and  the  total  stack  potential  are  monitored, 
the  produced  water  is  condensed  at  the  outlet  in  order  to  measure  its  amount.  The  model  used  is  a  mono-dimensional  description  of  a  single  cell. 
It  is  based  on  a  semi-empirical  law  for  the  local  electrochemical  response  computation.  The  local  pressures,  temperatures,  gas  compositions  and 
membrane  hydration  are  computed  dynamically.  A  semi-empirical  electrochemical  model  of  the  active  layer  is  fitted  to  experimental  data,  the  gas 
diffusion  layer  porosity  and  electro-osmosis  coefficient  of  the  membrane  have  been  adjusted  to  obtain  the  best  fit.  The  simulations  calculate  the  cell 
potential  and  the  current  density  profile  in  the  cell  for  various  feeding  conditions.  The  simulations  allow  us  to  explain  the  effect  of  gas  hydration 
on  the  stack  potential  and  to  predict  possible  drying  out  of  the  membrane. 
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1.  Introduction 

One  of  the  tasks  of  the  PICOS  project  is  to  characterize  the 
operation  of  a  1 .5  kW  Sorapec-manufactured  stack  of  12  cells  of 
330  cm2.  The  membranes  are  made  of  Nafion®  1 15,  we  have  no 
information  about  the  electrodes.  The  objective  of  this  charac¬ 
terization  is  the  evaluation  of  the  internal  temperatures,  partial 
pressures,  membrane  hydration,  liquid  water  saturation,  and 
resulting  current  density  distribution  during  operation  under  dif¬ 
ferent  experimental  conditions.  This  task  can  hardly  be  achieved 
by  direct  measurement,  because  adapted  sensors  (for  instance  to 
measure  the  partial  pressure  of  oxygen  inside  the  electrodes)  are 
not  available,  and  if  they  are,  it  is  not  possible  to  insert  them 
in  an  already  existing  and  closed  industrial  stack.  So  exper¬ 
imental  data  is  limited  to  fluids  inlet  and  outlet  conditions, 
global  stack  current  and  cells  potentials,  bipolar  plates  temper¬ 
atures.  The  way  to  go  ahead  of  these  practical  limitations  is  to 
build  a  model  of  what  happens  inside  the  cells.  To  be  able  to 
give  some  insight  on  the  internal  mechanisms,  this  model  has 
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to  be  based  on  physical  phenomena  involved  in  the  fuel  cell 
response. 

2.  Experimental  aspects 

The  stack  is  tested  on  a  specific  test  bench  (Fig.  1)  that 
allows  to  impose  independently  the  oxygen  (or  air)  and  hydro¬ 
gen  flow  rates,  the  gases’  temperature  and  water  content  at  the 
stack  inlet,  the  gas  pressure  at  the  stack  outlets  and  the  delivered 
current. 

The  stack  temperature  is  controlled  by  varying  the  tempera¬ 
ture  and  the  flow  rate  of  de-ionized  water  in  a  separated  cooling 
loop.  Thanks  to  a  fluid  distributor  inside  the  stack,  the  bipolar 
plates  are  all  cooled  in  parallel,  but  the  temperature  can  only 
be  regulated  on  one  plate  (the  one  located  in  the  middle  of  the 
stack  has  been  chosen).  Fig.  2  shows  that  the  temperatures  of 
the  other  plates  (excepted  end  plates)  are  close.  The  potential 
of  each  cell,  the  total  stack  potential  and  the  temperatures  of 
each  bipolar  plate  and  end  plates  are  monitored.  The  produced 
water  is  condensed  and  collected  at  the  outlet  and  its  amount  is 
measured  at  the  end  of  a  operation  cycle  in  order  to  estimate  the 
water  transfer  across  the  membrane  and  to  check  the  gas  humid¬ 
ification.  The  water  transfer  across  the  membrane  is  found  to 
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Fig.  1.  Stack  mounted  on  the  experimental  test  bench. 


be  negligible  compared  to  the  amount  of  water  produced  during 
operation. 

In  order  to  fit  the  model  (see  Section  4),  the  steady-state  stack 
voltage  is  recorded  for  different  values  of  temperature  and  pres¬ 
sure  (Fig.  3).  In  these  experimental  runs,  the  stack  is  fed  with 
pure  and  dry  oxygen  and  hydrogen  and  with  stoichiometric  ratios 
of  1 .4  for  oxygen  and  1 .05  for  hydrogen.  The  delivered  current  is 
successively  fixed  to  50,  100,  200  and  300  A.  Each  point  is  hold 
during  at  least  15  min  (depending  on  observed  deviation  of  cell 
voltage  in  mV  min-1)  in  order  to  be  sure  of  steady-state  opera¬ 
tion.  No  effect  of  gas  hydration  is  noticed  under  these  working 
conditions,  the  reproducibility  of  the  tests  is  good  and  all  the 
cells  have  roughly  the  same  voltage  (Fig.  2).  This  confirms  that 
there  is  no  fluid  misdistribution  in  the  stack  and  allows  (for 


Fig.  2.  Typical  voltage  and  temperature  distribution  in  the  stack  (H2/O2)  at 
200  A. 


Stack  temperature  (°C) 


Effect  of  pressure 


Fig.  3.  Influence  of  the  temperature  and  pressure  on  the  stack  voltage  for  dif¬ 
ferent  values  of  current. 

modeling  purpose)  to  represent  the  stack  by  an  equivalent  single 
cell. 

These  experimental  data  allow  to  fit  the  model,  mainly  the 
parameters  for  the  electrochemical  response  of  the  cell  (30). 
In  a  second  phase,  experimental  runs  using  air  instead  of  pure 
oxygen  allow  us  to  fit  the  porosity  of  the  gas  diffusion  layer  and 
to  observe  an  effect  of  gas  hydration  (Section  5). 

3.  Fuel  cell  model 

Our  model  results  of  several  PhD  and  project  studies  around 
polymer  electrolyte  fuel  cell  simulation  [1,3].  The  fuel  cell 
modeling  activity  at  CEA  started  about  10  years  ago  with  a  sta¬ 
tionary  model  able  to  simulate  the  steady-state  current-potential 
response  of  a  single  fuel  cell  [3].  This  model  based  on  [5]  was 
improved  in  order  to  take  into  account  two-phases  flows  in  the 
gas  diffusion  layers  and  in  the  active  layers  of  the  electrodes.  The 
main  purpose  of  this  work  was  to  explain  the  potential  drop  of 
the  cells  at  high  current  densities,  giving  a  qualitative  and  quan¬ 
titative  answer  to  the  question:  diffusion  limitation,  drying  out  of 
the  anode  or  flooding  of  the  cathode?  With  growing  interests  for 
transportation  applications  (implying  transient  operation)  and 
improvement  of  computing  techniques,  a  dynamic  model  has 
been  developed  [1]  in  order  to  simulate  non- stationary  cycles  [6], 
to  test  regulation  strategies  [7]  and  to  be  able  to  exploit  transient 
experimental  results  for  understanding  global  fuel  cell  opera¬ 
tion.  As  commercial  membrane-electrodes  assemblies  (MEA), 
without  access  to  internal  compositions  but  with  high  repeata¬ 
bility  are  now  used,  the  model  of  active  layer  is  simplified  to  a 
semi-empirical  law  fitted  on  experimental  data  [2] . 
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3.1.  General  description 

This  semi-empirical  law  allows  to  compute  the  potential  ver¬ 
sus  current  response  of  the  ME  A  given  the  local  pressures, 
temperatures  and  species  concentrations  at  the  active  layers 
[1,2].  These  local  conditions  are  dynamically  computed  using 
conservation  and  mass  transport  laws  in  the  different  parts  of  the 
cell  (channels,  gas  diffusion  layers,  membrane).  The  described 
phenomena  are:  pressure  drop  in  the  gas  distributors,  gas  diffu¬ 
sion  in  the  electrodes,  liquid  water  transport  in  the  electrodes, 
water  transport  in  the  membrane,  ohmic  loss  in  the  membrane, 
heat  production  and  transfer.  The  main  assumptions  supporting 
the  model  are  that  the  local  electrochemical  cell  response  time  is 
fast  (typically  less  than  that  10~3  s)  and  that  the  mass  and  ther¬ 
mal  transport  are  slow  (typically  from  10_1  to  102  s),  so  they 
have  to  be  computed  dynamically  in  the  context  of  transporta¬ 
tion  applications. 


3.2.  Gas  distributors 


Gas  distributors  (in  our  case  channels)  are  ID-meshed  from 
the  gas  inlets  to  the  gas  outlets.  Mass  and  thermal  balances  are 
computed  on  each  mesh  considering  gas  inlet,  gas  outlet,  oxygen 
or  hydrogen  consumption,  water  production,  water  condensation 
or  evaporation  and  heat  exchanges  with  the  bipolar  plate  and  the 
gas  diffusion  layer.  Pressure  drop  is  modeled  using  simplified 
linear  correlations  (1)  and  (2). 


Ggaz  — 
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£  ^  ^gaz 

0) 

AsCny 

^  APliq 

(2) 

These  relations  are  given  in  [8],  Q  is  the  volumic  flow  rate;  A  and 
L,  the  cross  section  and  the  length  of  the  channel,  respectively; 
s  the  fraction  of  volume  occupied  by  the  liquid  water;  Cq  a 
coefficient  to  fit  and  A P  the  pressure  drop.  The  liquid  phase  is 
represented  by  subscript  liq,  the  gas  phase  by  gaz. 


3.3.  Electrodes  ( hydraulic  aspect ) 


3.3.1.  Transport  phenomena 

The  diffusion  fluxes  in  the  gas  phase  are  driven  by  the  inverted 
Stefan-Maxwell  equations  (3)-(5).  The  molar  fluxes  F/  are 
computed  from  the  molar  fraction  difference  AX/,  taking  into 
account  the  thickness  e  and  the  surface  S  of  the  electrode,  the 
molar  masses  M/  of  the  species,  the  gas  concentration  Cg  and 
the  gas  density  pg  [9]. 
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The  diffusion  coefficients  D  Jj  in  a  ternary  (three  species)  mixture 
used  in  these  equations  are  deduced  from  the  binary  diffusion 
coefficients  Dy,  considering  the  molar  fraction  Xf  and  the  molar 
mass  Mi  of  each  species. 
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The  coefficients  Dy  used  in  the  porous  electrodes  are  corrected 
from  their  nominal  values  Dy  [1,3]  to  take  into  account  porosity 
£  of  the  material  and  liquid  water  saturation  s. 

DU  =  [(1  -  s)ep2Dij  (7) 

These  equations  are  completed  by  convective  transport  terms  for 
gas  and  liquid  phases  (8),  (9). 

h  =  KK rg(*)  A/Jg  (8) 

S  Mgvg  e 

F\  =  KKr\(s)  A  Pi 
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where  Krg(s)  =  (1  —  s)3  and  Xri(s)  =  s3  are  the  relative  permeabil¬ 
ities. 


The  hydraulic  model  of  electrode  is  decomposed  into  two 
parts,  the  one  is  computing  the  gas  and  water  transport  across 
the  gas  diffusion  layer  (GDL),  the  other  is  computing  the  mass 
and  heat  balance  inside  the  electrode.  A  simplified  bondgraph 
representation  [4]  of  this  approach  is  presented  in  Fig.  4. 
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Fig.  4.  Hydraulic  model  of  the  electrode. 


3.3.2.  Balance  equations 

The  model  computes  the  mass  balance  (10)— (13)  for  each 
component  (written  here  for  the  cathode  side): 


(10) 
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where  nspc  is  a  number  of  moles  of  the  specie  spc,  J]  Espc  the 
algebraic  sum  of  the  inlet  and  outlet  fluxes  for  this  specie,  Fy^\ 
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is  the  algebraic  rate  of  condensation  for  the  water  (negative  if 
water  evaporates). 

Vapor  saturation  is  then  calculated  (14),  (15)  and  rates  of 
evaporation  or  condensation  necessary  to  stay  close  to  the  liq¬ 
uid/vapor  equilibrium  are  estimated  (16),  (17). 
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else 

F\->1  =  —  ^  ^  Giq 

Gas  and  liquid  pressures  are  then  deduced,  considering  the  vol¬ 
ume  occupied  by  liquid  water  (18),  (19),  using  the  perfect  gas 
hypothesis  and  considering  capillary  effect  for  water  pressure 


calculation  (19)-(22). 
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where  J(s )  is  the  so  known  Leverett  function  [10],  Rv  the  pores 
sizes,  a  the  superficial  tension  and  0  contact  angle. 

At  last,  the  energy  balance  equation  is  used  to  compute  the 
temperature  in  each  mesh  (23): 
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Fig.  5.  Model  of  water  transport  in  the  membrane. 


mass  balance.  The  water  content  of  the  electrolyte  is  computed 
dynamically  from  the  interfacial  conditions. 

Equilibriums  with  hydration  conditions  in  the  electrodes  are 
given  [3,5]  for  anode  and  cathode  side  by  (24). 


if  a  <  1  then  0.043  +  17.81a  —  39.85a2  +  36.0a3 


else  14  +  1.4(a  —  1) 


(24) 


where  a  is  the  water  activity  in  the  electrode  at  the  interface  with 
the  membrane  and  X  is  the  water  content  of  the  membrane  at  the 
interface  with  the  electrode.  The  diffusion  coefficient  for  each 
mesh  is  then  estimated  [3]  using  (25) 


Dx  =  (6.707  x  10-8A.  +  6.387  x  10"7)exp 


(25) 


and  the  diffusion  flux  through  the  mesh  is  computed  using  (26) 
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where  T^diff  is  the  flux  of  diffusion  (mol  s-1)  and  em  is  the  thick¬ 
ness  of  the  half  membrane.  In  a  same  way,  the  electro-osmosis 
flux  is  given  [3]  by  (27) 
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where  Feo  is  the  molar  flux,  I  is  the  current  and  F  is  the  Faraday 
constant. 


3.5.  Electrochemical  response 


+  ^  F\igCpliq(Te  T ) 
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(23) 

where  TQ  is  the  inlet  temperature  for  each  component  and  T 
the  mixture  temperature  inside  the  electrode.  Cp  is  the  thermal 
capacity  and  Lv  the  latent  heat  for  water  evaporation.  The  contri¬ 
bution  V  ^  related  to  adiabatic  compression  is  usually  neglected 
because  we  do  not  study  very  sharp  pressure  transients. 


The  relations  giving  the  thermodynamic  potential  and  the 
activation  over-voltage  can  be  found  in  [1]  and  [2].  Global  cell 
potential  (for  one  mesh)  is  given  by  (28). 

U  —  £rev  T  77 act  FmI  (28) 

where 

£rev  =  ai  +  a2(T  -  298.15)  +  <*3r(0.5  lnPo2  +  lnPH2)  (29) 
is  the  thermodynamic  potential  and 

77 act  =  Pi  +  foT  +  P3T In  i  +  faT In  Pq2  (30) 


3.4.  Membrane 

The  model  of  the  membrane  considers  two  meshes  for  mass 
transport,  anode  and  cathode  sides  (Fig.  5),  and  one  mesh  for 


is  the  activation  over-voltage.  The  temperature  T  and  the  partial 
pressures  Pq2  ,  Pu2  are  the  local  conditions  at  the  active  layers  of 
the  electrodes.  The  coefficients  oti  are  calculated  from  physical 
parameters  [1,2].  The  coefficients  fr  are  fitted  to  experimental 
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results  (Section  4).  Rm  is  the  membrane  electrical  resistance.  The 
local  current  density  i  depends  on  the  amount  of  liquid  water  in 
the  electrode  according  to  (31) 


where  I  is  the  surface  current  density  and  s  the  volume  ratio  of 
liquid  water  in  the  electrode  (s  can  theoretically  vary  from  0  to  1 , 
simulations  have  given  typical  values  of  0.2  in  usual  conditions). 

At  last,  the  superficial  heat  production  is  given  [1]  by  (32). 

—  =  (1.2517  -  U)I  (32) 

s 

By  convention  in  this  model,  the  heat  production  is  calculated 
using  the  lower  heating  value  (LHV),  i.e.  considering  a  vapor 
water  production  (possible  condensation  that  depends  on  local 
temperature  and  pressure  conditions  is  considered  in  another  part 
of  the  model).  The  LHV  is  computed  for  the  standard  conditions 
(the  variation  with  the  temperature  or  with  the  pressure  is  lower 
than  1%  in  our  range  of  conditions  [1]). 

4.  Parameters  estimation 

The  parameters  that  have  to  be  identified  are  the  coefficients 
of  the  local  electrochemical  response.  To  be  able  to  identify 
these  parameters,  it  is  necessary  to  vary  current  density,  gases’ 
partial  pressure,  and  temperature.  One  difficulty  is  that  the  inde¬ 
pendent  variables  involved  in  the  model  are  different  from  those 
controlled  in  the  experiment:  the  experimental  setup  does  not 
allow  to  directly  impose  the  local  conditions,  but  they  have  to 
be  computed  starting  from  measured  parameters  (plate  tempera¬ 
ture,  pressure  in  the  gas  distributors,  inlet  gas  composition)  and 
using  the  model.  So  the  method  for  model  identification  is  itera¬ 
tive:  at  first,  we  compute  local  conditions  using  parameters  from 
the  literature  [1,2]  for  electrochemical  response.  Using  these 
computed  conditions  (hydrogen  and  oxygen  partial  pressures, 
local  current  density  and  membrane  resistance)  we  are  able  to 
fit  the  parameters  of  (30).  We  can  then  compute  again  the  local 
conditions,  using  the  updated  set  of  parameters.  These  steps  are 
iterated  until  the  parameters  do  not  change  any  more  from  one 
iteration  to  another.  In  fact,  the  computed  local  conditions  were 
found  not  to  be  very  sensitive  to  the  electrochemical  response, 
and  usually  two  iterations  were  enough  to  achieve  identification. 

The  parameters  fitting  is  performed  using  Matlab@’s  linear 
regression  capabilities.  Considering  (30)  we  define: 

xi  =  T,  x2  =  Tlog(/),  *3  =  riog(Po2)>  and 
X  =  [  1  x\  X2  X3]  (33) 

(where  x;  is  a  vector  in  which  number  of  elements  corresponds 
to  the  number  of  experiments  to  simulate)  to  obtain  a  linear 
relation  (34).  If  y  is  the  measured  overvoltage,  the  coefficients 
b  =  [  P2  P3  P4  ]  of  (30)  can  be  fitted  by  linear  regression 
of  the  model  (34). 

(34) 


Fig.  6.  Maximal  error  between  experiment  and  simulation  for  varying  electro¬ 
osmosis  coefficient. 


A  total  of  36  experimental  runs,  covering  the  most  working  con¬ 
ditions,  are  used  for  this  identification.  As  the  model  of  cell  is 
meshed  into  10  parts  from  inlet  to  outlet,  each  vector  jq  is  com¬ 
posed  of  360  values.  The  values  of  the  parameters  found  are 
(35): 

Pi  =  -8.25  x  10_1,  p2  =  1-87  x  10-3, 

P3  =  -2.16  x  10"4,  £4  =  8.75  x  10~5  (35) 

Another  iteration  level  was  necessary  to  find  out  the  optimal 
value  (the  one  that  allows  the  best  fit  for  the  other  parameters) 
of  the  electro-osmosis  coefficient  used  in  (27)  (Fig.  6). 

In  a  second  phase,  experimental  data  obtained  with  air  (with 
a  stoichiometric  ratio  fixed  to  2)  allowed  to  fit  the  porosity  of 
the  electrodes,  influencing  diffusion  coefficients  used  in  (3)-(5). 
The  parameters  of  the  electrochemical  response  were  of  course 
kept  unchanged.  A  porosity  of  0. 15  was  found  for  the  electrodes. 
This  value  is  low  compared  to  porosities  usually  found  in  liter¬ 
ature  [3,5],  but  we  can  note  that  the  porosity  is  related  to  the 
thickness  of  the  electrodes,  that  is  unknown. 
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Fig.  7.  Relative  humidity  in  the  channels  and  inside  the  electrodes  at  70  °C  and 
2.5  bar,  working  with  dry  H2  and  air.  Position  0  is  the  inlet  of  the  stack,  position 
1  is  the  outlet. 
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Relative  humidity 


Fig.  8.  Relative  humidity  in  the  channels  and  inside  the  electrodes  at  80  °C  and 
2.0  bar,  working  with  dry  H2  and  air.  Position  0  is  the  inlet  of  the  stack,  position 
1  is  the  outlet. 


Current  density  distribution 


Fig.  9.  Current  distributions  in  the  stack  for  different  conditions  (dry  gases). 
Position  0  is  to  the  inlet  of  the  stack,  position  1  is  the  outlet. 

5.  Stack  operation  characterization 

As  the  model  is  fitted,  we  can  perform  simulations  in  order 
to  access  local  gas  composition  and  current  distribution.  The 
simulations  show  that  for  the  stack  working  with  air  at  2.5  bar 
absolute  pressure,  the  membrane  will  not  dry  out.  This  can  be 
explained  by  Fig.  7,  the  vapor  diffusion  through  the  gas  diffusion 
layer  limits  the  water  transfer.  Consequently,  the  current  density 


I  (A.cnr2) 


Verticalposition  (mm) 

Fig.  10.  Experimental  current  distribution  for  H2/O2  gas  feeding.  Two  experi¬ 
mental  runs  are  compared. 


is  homogeneous  along  the  surface  of  the  cells,  and  no  effect  on 
the  cell  voltage  is  expected. 

If  we  decrease  the  working  pressure  down  to  2.0  while 
increasing  the  temperature  up  to  80  °C,  simulations  predict  a 
drying  out  of  the  membrane  and  active  layers  at  the  stack  inlet 
as  shown  in  Fig.  8.  This  leads  to  inhomogeneous  current  distri¬ 
butions  along  the  cell,  as  illustrated  in  Fig.  9. 

These  simulations  are  confirmed  by  experimental  results:  an 
effect  of  gas  hydration  on  stack  voltage  was  in  fact  observed 
at  2.0  bars  and  80  °C,  whereas  no  effect  was  noticed  at 
2.5  bars. 

6.  Conclusion 

We  have  developed  a  model  that  allows  to  access  inter¬ 
nal  operation  of  the  stack.  This  model  was  successfully  fitted 
on  experimental  data  obtained  on  a  test  bench  in  the  context 
of  the  PICOS  project.  The  fitted  model  is  able  to  compute 
the  one-dimensional  current  density  profile  in  the  stack  and 
to  reproduce  effects  of  temperature,  pressure  and  gas  humid¬ 
ification  on  the  stack  voltage.  Effects  of  low  gas  humidifica¬ 
tion  and  resulting  drying  out  of  the  membrane  (depending  on 
stoichiometric  ratios,  pressure  and  temperature)  can  be  pre¬ 
dicted.  Experimental  devices  allowing  to  perform  a  measure¬ 
ment  of  current  densities  (based  on  magnetic  field)  are  currently 
under  development  in  order  to  confirm  the  predictions  of  the 
model. 

First  tests  in  O2/H2  gas  alimentation  are  promising  (Fig.  10), 
the  measured  current  density  is  in  accordance  with  the  total  cur¬ 
rent.  As  expected  under  these  conditions,  the  current  distribution 
obtained  is  nearly  uniform.  These  results  can  be  compared  to 
simulation  results  (Fig.  9).  The  differences  may  be  attributed  to 
mechanical  constraints  in  the  stack  and  border  effects  that  are 
so  far  not  taken  into  account  in  the  model. 
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